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Abstract 

Oxidative stress induced by excessive production of reactive oxygen species (ROS) has been implicated in the etiology of 
many human diseases. It has been reported that fullerenes and some of their derivatives-carboxyfullerenes-exhibits a 
strong free radical scavenging capacity. The permeation of C 50 -fullerene and its amphiphilic derivatives-C 3 -tris-malonic-C 50 - 
fullerene (C 3 ) and D 3 -tris-malonyl-C 60 -fullerene (D 3 )-through a lipid bilayer mimicking the eukaryotic cell membrane was 
studied using molecular dynamics (MD) simulations. The free energy profiles along the normal to the bilayer composed of 
1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) for C 60 , C 3 and D 3 were calculated. We found that C 60 molecules alone 
or in clusters spontaneously translocate to the hydrophobic core of the membrane and stay inside the bilayer during the 
whole period of simulation time. The incorporation of cluster of fullerenes inside the bilayer changes properties of the 
bilayer and leads to its deformation. In simulations of the tris-malonic fullerenes we discovered that both isomers, C 3 and D 3 , 
adsorb at the surface of the bilayer but only C 3 tends to be buried in the area of the lipid headgroups forming hydrophobic 
contacts with the lipid tails. We hypothesize that such position has implications for ROS scavenging mechanism in the 
specific cell compartments. 
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Introduction 

Fullerenes and their derivatives are considered perspective for 
various applications in medicine and pharmacology [1]. Particu- 
larly, C 60 -fullerene and tris-malonic C 60 -fullerene (Figure 1) are 
known as free radical scavengers, that have been shown to protect 
cells from free radicals (including ROS) that can induce apoptotic 
injuries in vitro [2,3] as well as in different cell types: neuronal cells 
[4], hepatoma cells [5] and epithelial cells [6]. At the same time, 
some of fullerene derivatives (for example dendrofullerenes) [7] 
inhibit HIV-protease, what gives another prospect for their 
biomedical use. Harhaji et al. demonstrated [8] the antitumor 
effect of the water suspension of Cen on glioma cell cultures when 
irradiated by light: high concentrations of fullerene caused 
necrosis, while low concentrations stopped proliferation of the 
cells and eventually lead to autophagy. 

Fumelli et al. [9] have demonstrated that carboxyfullerenes 
protected human keratinocytes from apoptosis induced by 
ultraviolet-B (UVB) irradiation. Later the possibilities of prevent- 
ing the neurotoxicity generated by levodopa (L-3, 4-dihydroxy- 
phenylalanine) have been demonstrated as chemical properties of 
a water-soluble fullerene derivative and ascorbic acid [10]. Wang 
et al. compared [1 1] antioxidant activity of the two regioisomeric 
forms of tris-malonic fullerene (Figure 1)-C 3 (with the malonic 
groups localized on one side of the fullerene molecule) and D 3 
(with the malonic groups distributed evenly). It was found that C 3 
has a higher activity against ROS than D 3 , presumably due to its 



better interaction with biomembranes. However, the molecular 
basis of these interactions remains elusive. Particularly, the process 
of passive transport through cell membranes of Cgo and its 
derivatives must be elucidated. 

Most of the fullerenes get inside the cell by the energy- 
dependent mechanism of endocytosis [12]. On the other hand, 
hydrophobic C 6(l molecules can directly interact with the 
biomembrane and penetrate it in an energy-independent manner. 
One of the theoretical methods, which allows to tackle this 
problem, is computational simulation by molecular dynamics. 
These can be used to study the penetration of small molecules [13] 
as well as nanoparticles [14] through biomembranes. Several 
atomistic molecular dynamics simulation studies have been 
performed that investigate the interactions of C 60 and its 
derivatives with model biological membranes [15,16,17]. It was 
shown that water-soluble fullerene derivatives do not penetrate 
into the membranes as the C 60 does [16,17,18]. In [15] a coarse- 
grain model was applied to study penetration abilities of multiple 
fullerenes into different membranes and the authors showed that 
multiple fullerenes inside the bilayer cause no deformation of the 
membrane. 

Here we report a comparative study of the interaction of C 60 - 
fullerene and its amphiphilic derivatives, C 3 and D 3 , with 
eukaryotic membranes composed of DPPC lipids using molecular 
dynamics simulations. We performed series of unrestrained MD 
simulations and calculated free energy profiles as a function of the 
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Figure 1. The studied species of the C 60 fullerene. A. 

Buckminsterfullerene C 60 . B. D 3 stereoisomer of tris-malonyl-C 60 - 
fullerene. C. C 3 stereoisomer of tris-malonyl-C 60 -fullerene. 
doi:10.1371/journal.pone.0102487.g001 

position of the studied molecules along the bilayer normal. The 
aim of our investigation is to explain the mechanism of ROS 
scavenging for fullerene species with different activities. We found 
out that molecules have different penetration depths into the 
bilayer what has implications regarding the possible antioxidant 
effects of G 60 carboxyl derivatives. 

Materials and Methods 

Four different systems were used for the simulations: a lipid 
bilayer with one molecule C 60 , with ten molecules of C 60 , with one 
C3 molecule and with one D 3 molecule. The model of the 
eukaryotic membrane consisted of 128 solvated DPPC lipids. The 
Ceo - , C3- or D 3 -molecules were placed in water and did not 
contact the lipid surface in the beginning of the simulations. In all 
simulations we used the Berger model for lipids [19] and the SPC/ 
E water model [20]. The force field parameters for C 60 , C 3 and D 3 
were set as follows. The lengths of the two types of C-C bonds in 
fullerenes were set to 0,139 nm and 0,144 nm according to [21]. 
Parameters of the van der Waals interactions were taken from 
[16,17]. The atomic partial charges of carbon atoms in a fullerene 
were set 0 for the all of them due to the symmetry of the molecule, 
while the Mulliken partial atomic charges for C 3 and D 3 
derivatives were calculated in GAMESS suite (see SI, Table SI) 
using ab initio Hartree-Fock method with 6-3 1G** basis set 
(prefaced by the geometry optimization of the molecules at the 
same level of theory). 

The all of MD simulations were carried out using Gromacs 4.5 
[22]. Simulations were carried out in the NPT ensemble: semi- 
isotropic pressure coupling (Parrinnelo-Rahman barostat [23], 
time constant 2 ps) and constant temperature 323 K (Nose-Hoover 
thermostat [24,25], time constant 2 ps). Lipids and water were 
coupled independently to the heat bath. Periodic boundary 
conditions were applied in all three dimensions. All bond lengths 
were kept constant using the LINCS [26] algorithm. The time step 
was 2 fs, long-range electrostatic interactions were treated with the 
PME algorithm [27] (real space cutoff 1 nm, FFT grid spacing 
0.18 nm). The Lennard-Jones potentials were computed by using 
a cutoff length of 1.2 nm. 

Free energy profiles for fullerenes were obtained using the 
metadynamics approach [28] and the plugin PLUMED [29]. 
Each of the metadynamics calculations was preceded with an 
equilibration run. A total of 1.5 |is of molecular dynamics was 
carried out. We chose a z-component (perpendicular to the 
membrane plane) of a vector connecting centers of masses of the 
membrane and the fullerene as a reaction coordinate (called 
collective variable) for the metadynamics simulations. The following 
parameters were used for metadynamics: the integration time step 



was decreased to 1 fs to improve stability of the systems, the time 
interval between the addition of two Gaussian functions was 
t = 500 fs, the Gaussian height w = 0.5 kj/mol, and the Gaussian 
width d = 0.5 A. Potential walls were applied to keep the values of 
the chosen collective variables in the area of a single membrane 
leaflet with the adjacent water layer, i.e. between — 10 and 35 A. 
Simulations were carried on until the convergence of the free 
energy surfaces was achieved after 45 ns. Additional 4 indepen- 
dent trajectories were obtained for each of the studied systems to 
calculate the average free energy profiles along with standard 
deviation. Sufficient sampling of the collective variables in 
metadynamics simulations was checked by plotting their progress 
along the simulation time (Figure SI) to confirm adequate 
sampling of the chosen reaction coordinate. All of the performed 
MD runs are listed in Table S2. Analysis of the obtained 
trajectories was done with standard Gromacs 4.5 utilities. 
Membrane properties calculations were made using GridMAT- 
MD [30] . Plots of molecular hydrophobic potential (MHP) for the 
membrane surfaces were built using PLATINUM software [31]. 

Results and Discussion 

C 60 fullerene rapidly penetrates into membrane 

The equilibrium MD simulation of a single Ceo performed in 
our study showed that during the first nanosecond of the MD 
trajectory the fullerene was adsorbed in the area of the lipid 
headgroups. However, in the third nanosecond a spontaneous dive 
into the tail region was observed (Figure 2A, C), which is consistent 
with [15,16,17]. Over the following almost 100 ns of the 
simulation the fullerene remains inside the membrane. 

Table 1 represents a summary of previous studies of interactions 
of C 6 o and different membranes. According to this table we can 
conclude that the shape of the PMF profiles for the fullerene 
interacting with the membrane are quite similar. All methods, 
which were used to describe interactions between C 60 and 
membranes, show a low minimum inside the membrane (near 
its center). Yet the details are different: some simulations show a 
small barrier for fullerene before entering the tail region. Also it's 
important to notice that coarse-grain (CG) models underestimate 
free energy values for Cgo Different methods of free energy 
calculations compromise energies: due to the small simulation time 
it is almost impossible to reach the thermodynamic equilibrium in 
such inhomogeneous systems as membranes. 

We have used the metadynamics approach to calculate the free 
energy profiles of the fullerenes interacting with the biomembrane, 
which, to the best of our knowledge, has not been applied to the 
fullerene-membrane systems before. So we have validated the 
method against the well-studied case of the C 60 fullerene 
interacting with the membrane and came up with the conclusion 
that the metadynamics can be efficiently exploited in case of the 
C 3 and D 3 derivatives as well. In the figure 2B one can see that the 
fullerene overcomes the energy barrier of approximately 15 kj/ 
mol to pass into the hydrophobic region of the lipid tails. The 
global energy minimum for the G 60 is located 0.7 nm behind the 
center of the membrane: if the Cgo is moving closer to the center of 
the membrane its free energy is increasing by almost 30 kj/ mol. 
So, the Ggo remains at a certain distance from the precise 
membrane center in the region of the hydrocarbon chains. To 
double-check this observation Cgo was placed in the center of the 
membrane (between C-terminal atoms of lipid tails) and a 50 ns 
simulation was run. The fullerene in this simulation also migrated 
to the area located approximately 0,8 nm from the center of the 
membrane and remained there for the rest of the simulation. 
Absorption of a single C so on the surface of membrane and its 
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Figure 2. Characteristics of fullerene-membrane interactions. A. Distance between the center of mass (COM) of C 60 and COM of the 

membrane. On the third nanosecond fullerene spontaneously "jump" into the membrane (the membrane surface is shown with the dot line). B. Free 
energy profile of the process of the C 60 penetration into the model eukaryotic membrane. Potential wall at 30 A is shown as the dotted line. C. 
Snapshot of the system with a single C 60 molecule inside the membrane. 
doi:1 0.1 371 /journal.pone.01 02487.g002 



"jump" to the tails region cause no damage to the membrane (see 
Table 2). The above-mentioned results appear to be in a good 
agreement with the previously reported data serving as a good 
validation of the chosen model. 

Cluster of fullerenes causes changes in the physical 
properties of the membrane upon penetration 

In an aqueous medium fullerenes exist only in aggregates 
because of their hydrophobic nature. In our study we provide 
calculations of ten fullerenes and a DPPC membrane with heavy- 
atom resolution. To study the permeability of the membrane for a 
fullerene cluster which exists in real solutions, we assembled a 
system consisting of 10 fullerenes above the surface of the 
membrane. At the start of the simulation fullerenes did not 
contact with each other or with the membrane, but after a few 
nanoseconds, they aggregated and absorbed in the area of the lipid 
headgroups. The first fullerene from the cluster "jumped" inside 
the membrane after 3 ns. The rest of the fullerenes penetrated the 
membrane afterwards, so that eventually, after 100 ns, 9 of 10 
fullerenes had dived into the membrane (Figure 3C). All of them 
stayed in the membrane during the next 400 ns of the simulation 
time. We did not observe their disaggregation during simulation 
time but according to [15] this happens in the microsecond scale. 
The most probable localization of fullerenes inside the membrane 
is around 0.5-1 nm from the center. Contrarily to the previously 
reported coarse-grained MD simulations [15], the membrane is 
deformed significandy by the addition of a large number of 
fullerenes: it curves (Figure 3C, average curvature radius 
calculated as in [32], Figure S2), the area per lipid head reduces 
from 69 to 56 A 2 , on the contrary, thickness increases to 38,7 A 
and up to 42 A at the end of simulation. Such changes are 
connected to the formation of a fullerene cluster inside the 
membrane (Figure 3C) consisting of 8 fullerenes. MHP plots were 



made for both membrane sides of the two membrane conforma- 
tions: at the beginning and at the end of the simulation (see SI: 
Figure S3). We observed that a large hydrophobic area had 
appeared on the membrane surface from the side where the 
fullerenes penetrated the membrane (Figure 3D). It means that 
when the fullerenes get inside the membrane a small pore appears 
with the hydrophobic lipid tails being exposed to the solvent. 
Lipids in this conformation are more vulnerable to ROS. The 
diffusion coefficients for DPPC molecules do not change upon 
addition of the fullerenes in the membrane and they are in a good 
agreement with the experimental data [33]. Pore or micelle 
formation is not observed though possibly due to insufficient 
simulation time. The order parameters for acyl chains Sqh 
(Figure 3 A and 3B) for lipid chains sn-1 and sn-2 were calculated 
using the following equation: 



S c „=--(3cos 2 e CH -\) 



From these plots one can see that the order parameters do not 
change significantly when only one fullerene is inside the 
membrane while they modestly increase in the simulation with 9 
fullerenes inside the membrane. 

Our observations indicate that the addition of multiple C 60 
molecules does not cause a phase transition of the lipid bilayer but 
leads to its loosening and bending as well as to the formation of 
hydrophobic pores, what altogether can increase membrane 
susceptibility to ROS and explains toxicity of high concentrations 
of the C6o, which was reported previously [34]. 
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Table 2. Properties of the DPPC membrane in simulations with different amount of fullerenes. 





Membrane property 


Pure DPPC 


DPPC with one fullerene 


DPPC with ten fullerenes 










70-100 ns of simulation 


470-500 ns of simulation 


Thickness, A 


33,7 


32,0 


38,7 


42,2 


Area per lipid, A 2 


69 


71 


56 


56,6 


Lateral diffusion coefficient, 10~ 7 cm 2 /s 


0,98 


1,3 


0,99 


0,61 


Membrane curvature 






32 


30,2 



doi:1 0.1 371 /journal.pone.01 02487.t002 



C 3 and D 3 fullerene derivatives do not penetrate 
membrane 

We have carried out equilibrium MD simulations of the C3 and 
D 3 fullerenes interacting with the membrane along with 



metadynamics simulations in order to obtain a complete picture 
of molecular and energetical details of the membrane insertion 
processes for the C3 and D 3 variants (Figure 4, 5). 

The fully deprotonated forms of both studied fullerene 
derivatives should be the most probable species in aqueous 
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Figure 3. Order parameters, a snapshot of the final conformation of the system with ten C 50 molecules and MHP map of the final 
conformation. A-B. Order parameters for lipid tails calculated from the equilibrium MD trajectories for the pure DPPC membrane, DPPC with one 
C 60 molecule and DPPC with 10 C 60 molecules: A-sn1-chain, B-sn-2 chain. C. Snapshot of the final conformation of the system with ten C 60 molecules 
(totally nine fullerenes penetrated into the membrane). D. Molecular hydrophobic potential (MHP) map made for the upper surface of DPPC 
membrane. Carbone atoms of C 60 are rendered with red VdW spheres. 
doi:1 0.1 371 /journal.pone.01 02487.g003 
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C 3 - membrane COM distance, nm 
3.5 3.0 2.5 2.0 1.5 1.0 0.5 0 




Figure 4. C 3 stereoisomer of tris-malonic fullerene. A. Free energy profile of the process of the C 3 penetration into the model eukaryotic 
membrane. B. Intermediate orientation of the C 3 molecule adsorbed to the membrane with its solvent shell retained. C. Stable conformation 
(corresponding to the global energy minimum of the free energy profile) of C 3 adsorbed to the membrane and established hydrophobic contact with 
the lipid tails region. 
doi:10.1371/journal.pone.0102487.g004 



environment (pK a [=2.3 and pK a 2 = 5.7 for malonic acid in 
water). This seems to be no longer the case in the membrane 
hydrophobic region, where pK a s of acidic groups shift towards 
higher values. Higher pK a s values should favor protonated species 
and reduce thus the free energy cost for the translocation of the 
C 3 /D 3 to the region of the lipid tails. On the other hand in the 
region of the zwitterionic lipid head groups pK a of malonyl groups 
depends on the interplay of the local interactions of the fullerene 
derivatives with phosphate or choline groups. However, we 
assumed that these effects play secondary role at the stage of the 
C 3 /D 3 absorption to the bilayer and performed simulations of the 
fully ionized C 3 /D 3 species (charged —6) only. 

In the equilibrium simulation with C 3 , the C 3 molecule 
adsorbed to the membrane surface in two steps: in the beginning, 



it was localized at a distance of 2 nm from the COM of the 
membrane and retained its solvent shell (which consisted of 40 
water molecules) as seen in Figure 4 and, later, it moved closer to 
the COM of the membrane by approximately 0.5 nm and 
established a hydrophobic contact with the lipid hydrocarbon tails 
as well as ionic bonds with the lipid headgroups (Figure 4). This 
notable behavior of the C 3 derivative takes root in its 
stereochemistry-the molecule has two hemispheres: one is 
hydrophobic and the opposite carries three malonic acid residues 
and thus is hydrophilic (Figure 1C). We calculated the orientation 
of C 3 molecule in respect to the membrane normal along the 
equilibrium trajectory of the C 3 absorption to the membrane 
surface. Since the absorption process was completed after 25 ns, 
the C 3 upheld a specific orientation with its hydrophobic 



D 3 - membrane COM distance, nm 

3.5 3.0 2.5 2.0 1.5 1.0 0.5 0 




-160- 



Figure 5. D 3 stereoisomer of tris-malonic fullereneB A. Free energy profile of the process of the D 3 penetration into the model eukaryotic 
membrane. B. Orientation (corresponding to the global energy minimum of the free energy profile) of the D 3 molecule adsorbed to the membrane. 
doi:1 0.1 371 /journal.pone.01 02487.g005 
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hemisphere pointing to the middle of the membrane and the 
malonic groups interacting with the head groups of the lipids 
remained accessible to the solvent (Figure S4). Energetically, this 
orientation of the C 3 corresponds to the global minimum of the 
free energy profile. 

For both the malonic derivatives of Cgo we also calculated their 
solvent accessible surface area (SAS). Upon the absorption of C 3 to 
the membrane the SAS of C 3 decreased significantly (Figure S5). 
On the other hand, the SAS of D 3 remained almost the same 
during the whole simulation time, what indicates the lack of a close 
hydrophobic contact with the lipid tails in this case. 

The observed orientation of C 3 upon its absorption to the 
membrane surface appears to play a crucial role in efficient 
scavenging of ROS and radical protection: It resides not inside the 
membrane but near the lipid tail region, preventing chain 
peroxidation of the lipids. The isomeric form of C 3 -D 3 - does 
not have [35] such high cytoprotective properties since its three 
malonic acid residues cover the fullerene surface symmetrically 
(Figure IB) and D 3 only adsorbs efficiently to the region of the 
lipid headgroups where the distance from the COM of the 
membrane equals 2 nm and retains its solvent shell (see Figure 5C). 

Neither the C 3 nor the D 3 stereoisomers penetrate the 
membrane because of the large free energy penalty of being in 
the middle of the membrane for these molecules (Figure 4A and 
5A). 

Conclusions 

The processes of penetration and accumulation of C 60 fullerene, 
its cluster and the two derivatives of C 60 in lipid bilayers were 
studied using equilibrium MD and metadynamics approach. 
Absorption of all of the studied fullerene species on the DPPC 
membrane occurs within the first nanoseconds. C 60 fullerene (a 
single molecule and a cluster of ten molecules) spontaneously 
penetrates the membrane and remains inside of it throughout the 
whole simulation time. After penetration into the eukaryotic 
membrane, Cgo stays at a distance of 0.7— 0.8 nm from the center 
of the membrane. A cluster of fullerenes deforms the membrane 
causing its curvature and formation of hydrophobic regions, which 
presumably make the membrane more susceptible to the ROS 
attacks. 

The studied tris-malonic derivatives of the C 60 fullerene do not 
penetrate membrane but rather are accumulated at its surface. 
However, during adsorption to the membrane the C 3 , in contrast 
to the D 3 , appears to be deeply buried in the area of the lipid 
headgroups of the membrane forming specific hydrophobic 

References 

1 . Bosi S, Da Ros T, Spalluto G, Prato M (2003) Fullerene derivatives: an attractive 
tool for biological applications. Eur J Med Ghem 38: 913-923. 

2. Lin AMY, Chyi BY, Wang SD, Yu HH, Kanakamma PP, ct al. (1999) 
Carboxyfullercne prevents iron-induced oxidative stress in rat brain. Journal of 
Neurochemistry 72: 1634-1640. 

3. Lin AMY, Fang SF, Lin SZ, Chou CK, Luh 'FY, ct al. (2002) Local 
carboxyfullercne protects cortical infarction in rat brain. Neuroscicnec Research 
43: 317-321. 

4. Dugan LL, Turetsky DM, Du C, Lobncr D, Wheeler M, ct al. (1997) 
Carboxyfullerenes as neuroprotective agents. Proceedings of the National 
Academy of Sciences of the United States of America 94: 9434—9439. 

5. Huang YL, Shcn CKF, Luh 'FY, Yang HC, Hwang KG, ct al. (1998) Blockage 
of apoptotie signaling of transforming growth factor-beta in human hepatoma 
cells by carboxyfullercne. European Journal of Biochemistry 254: 38—43. 

6. Straface E, Natalini B, Monti D, Franceschi G, Schcttini G, ct al. (1999) C 3 - 
fullero-tris-methanodicarboxylic acid protects epithelial cells from radiation- 
induced anoikia by influencing cell adhesion ability. 1 ; EBS Lett 454: 335-340. 

7. Marchcsan S, Da Ros T, Spalluto G, Balzarini J, Prato M (2005) Anti-HIV 
properties of cationic fullerene derivatives. Bioorg Med Chem Lett 15: 3615- 
3618. 



contacts with the hydrocarbon tails of the lipids along with the 
ionic interactions with the head region and hydrophilic contacts 
with the solvent. This energetically favorable orientation of the C 3 
derivative, as we assume, accounts for its high ROS scavenging 
activity comparing to the D 3 . 

Supporting Information 

Figure SI Progress of the collective variable during the 
^605 C 3 and D 3 metadynamics simulations. 

(TIF) 

Figure S2 Membrane curvature in xz and yz directions. 
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